Change the positions of the galaxies in a SHAM to be shuffled then NFW-distributed, instead of on the subhalos.

Shuffling procedure is as followed, from Jeremey

procedure:

take a bin in halo mass (small bins, like 0.1dex wide). (this is all halos, regardless of whether they have a galaxy in them or not). take all the centrals and put them in a list. take all the satellites and put them in a separate list.

randomly assign the centrals to all the halos in the bin.

randomly assign each satellite to a halo in the bin (repeat until all satellites are gone. this should preserve poisson distribution of satellite occupation). when assigning a satellite to a halo, preserve the position of the satellite and velocity of the satellite relative to the original host halo. ie, your list of satllites has dx, dy, dz, and dvx, dvy, dvz, then you add x, y, z, and vx, vy, vz of the new halo to those quantities.



In [1]:

    
import numpy as np
import astropy
from itertools import izip
from pearce.mocks import compute_prim_haloprop_bins, cat_dict
from pearce.mocks.customHODModels import *
from halotools.utils.table_utils import compute_conditional_percentiles
from halotools.utils import *



In [2]:

    
from matplotlib import pyplot as plt
%matplotlib inline
import seaborn as sns
sns.set()



In [3]:

    
PMASS = 591421440.0000001 #chinchilla 400/ 2048
Lbox = 400.0
#catalog = np.loadtxt('ab_sham_hod_data_cut.npy')
catalog = astropy.table.Table.read('abmatched_halos.hdf5', format = 'hdf5')









    



WARNING: path= was not specified but multiple tables are present, reading in first available table (path=abmatched_halos.hdf5) [astropy.io.misc.hdf5]
WARNING:astropy:path= was not specified but multiple tables are present, reading in first available table (path=abmatched_halos.hdf5)



In [4]:

    
catalog = catalog[catalog['halo_mvir'] > 200*PMASS]



In [5]:

    
catalog.colnames









    Out[5]:





['halo_upid',
 'halo_y',
 'halo_x',
 'halo_z',
 'halo_rvir',
 'halo_vpeak',
 'halo_rs_klypin',
 'halo_snapnum',
 'halo_halfmass_scale',
 'halo_id',
 'halo_vx',
 'halo_vy',
 'halo_vz',
 'halo_rs',
 'halo_mvir',
 'halo_nfw_conc',
 'halo_vpeak_mag',
 'halo_vvir_mag',
 'halo_alpha_05_mag',
 'halo_shuffled_vpeak_mag',
 'halo_shuffled_vvir_mag',
 'halo_shuffled_alpha_05_mag',
 'host_halo_nfw_conc',
 'host_halo_rvir',
 'halo_nfw_x',
 'halo_nfw_y',
 'halo_nfw_z',
 'halo_sh_shuffled_vpeak_mag',
 'halo_sh_shuffled_vvir_mag',
 'halo_sh_shuffled_alpha_05_mag',
 'halo_shuffled_x',
 'halo_shuffled_y',
 'halo_shuffled_z',
 'halo_shuffled_upid',
 'halo_shuffled_host_mvir',
 'halo_hostid',
 'halo_x_host_halo',
 'halo_y_host_halo',
 'halo_z_host_halo',
 'halo_nfw_conc_host_halo',
 'halo_mvir_host_halo',
 'halo_rvir_host_halo',
 'halo_sh_shuffled_x',
 'halo_sh_shuffled_y',
 'halo_sh_shuffled_z',
 'halo_sh_shuffled_upid',
 'halo_sh_shuffled_host_mvir']

#del catalog['halo_shuffled_x'] #del catalog['halo_shuffled_y'] #del catalog['halo_shuffled_z'] #del catalog['halo_shuffled_nfw_conc'] #del catalog['halo_shuffled_rvir'] #del catalog['halo_shuffled_mvir'] del catalog['halo_shuffled_host_mvir'] del catalog['halo_shuffled_upid']

#del catalog['halo_sh_shuffled_x'] #del catalog['halo_sh_shuffled_y'] #del catalog['halo_sh_shuffled_z'] #del catalog['halo_sh_shuffled_nfw_conc'] #del catalog['halo_sh_shuffled_rvir'] #del catalog['halo_sh_shuffled_mvir'] del catalog['halo_sh_shuffled_host_mvir'] del catalog['halo_sh_shuffled_upid']



In [6]:

    
add_halo_hostid(catalog, delete_possibly_existing_column=True)



In [7]:

    
for prop in ['halo_x', 'halo_y', 'halo_z', 'halo_nfw_conc', 'halo_mvir', 'halo_rvir']:
    broadcast_host_halo_property(catalog, prop, delete_possibly_existing_column=True)



In [8]:

    
from halotools.utils.table_utils import compute_prim_haloprop_bins
from math import ceil
min_log_mass = np.log10(np.min(catalog['halo_mvir']))-0.001
max_log_mass = np.log10(np.max(catalog['halo_mvir']))+0.001
dlog10_prim_haloprop = 0.1
num_prim_haloprop_bins = (max_log_mass - min_log_mass) / dlog10_prim_haloprop
prim_haloprop_bin_boundaries = np.logspace(min_log_mass, max_log_mass,
    num=int(ceil(num_prim_haloprop_bins)))

prim_haloprop_bins = compute_prim_haloprop_bins(prim_haloprop = catalog['halo_mvir_host_halo'],\
                                                dlog10_prim_haloprop=dlog10_prim_haloprop,
                                                prim_haloprop_bin_boundaries = prim_haloprop_bin_boundaries)



In [9]:

    
shuffled_pos = np.zeros((len(catalog), 3))
shuffled_upids = np.zeros((len(catalog)))
shuffled_host_mvir = np.zeros((len(catalog)))



In [10]:

    
shuffled_mags = np.zeros((len(catalog), 3))
#shuffled_mags[:, 0] = catalog['halo_vpeak_mag']
#shuffled_mags[:, 1] = catalog['halo_vvir_mag']
#shuffled_mags[:, 2] = catalog['halo_alpha_05_mag']



In [11]:

    
from pearce.mocks import cat_dict
cosmo_params = {'simname':'chinchilla', 'Lbox':400.0, 'scale_factors':[0.658, 1.0]}
cat = cat_dict[cosmo_params['simname']](**cosmo_params)#construct the specified catalog!
cat.load_model(1.0, HOD = 'redMagic')



In [12]:

    
np.random.seed(64)
bins_in_halocat = set(prim_haloprop_bins)

for ibin in bins_in_halocat:
    
    if ibin==0:
        continue
    indices_of_prim_haloprop_bin = np.where(prim_haloprop_bins == ibin)[0]
    
    centrals_idx = np.where(catalog[indices_of_prim_haloprop_bin]['halo_upid'] == -1)[0]
    n_centrals = len(centrals_idx)
    satellites_idx = np.where(catalog[indices_of_prim_haloprop_bin]['halo_upid']!=-1)[0]
    n_satellites = len(satellites_idx)
    
    if centrals_idx.shape[0]!=0:
        rand_central_idxs = np.random.choice(indices_of_prim_haloprop_bin[centrals_idx], size = n_centrals, replace = False)
    else:
        rand_central_idxs = np.array([])

    for idx, coord in enumerate(['vpeak', 'vvir', 'alpha_05']):
        shuffled_mags[indices_of_prim_haloprop_bin[centrals_idx], idx]= \
                catalog[rand_central_idxs]['halo_'+coord+'_mag']
            
        shuffled_mags[indices_of_prim_haloprop_bin[satellites_idx],idx ] = \
                catalog[indices_of_prim_haloprop_bin[satellites_idx]]['halo_'+coord+'_mag']
    #Create second rand_central_idxs, Iterate through satellite hosts and assign them when they match. 
            
    #after randomly shuffling the satellites, we have to choose them a new host
    if centrals_idx.shape[0]!=0:
        rand_host_idxs = np.random.choice(indices_of_prim_haloprop_bin[centrals_idx], size = n_satellites, replace = True)
        #rand_host_idxs_sh = 
    else:
        rand_host_idxs = np.array([])
        #rand_host_idxs_sh = np.array([])
        
        
    hc_x, hc_y, hc_z = cat.model.model_dictionary['satellites_profile'].mc_halo_centric_pos(\
                                                        catalog[rand_host_idxs]['halo_nfw_conc'],
                                                        halo_radius = catalog[rand_host_idxs]['halo_rvir'])
    #hc_x = hc_y = hc_z = np.zeros_like(rand_host_idxs)
                 
    for idx, (coord, hc) in enumerate(izip(['x','y','z'], [hc_x, hc_y, hc_z])):
        shuffled_pos[indices_of_prim_haloprop_bin[centrals_idx], idx] = \
                catalog[indices_of_prim_haloprop_bin[centrals_idx]]['halo_'+coord]

                                        
        #shuffled_pos[indices_of_prim_haloprop_bin[satellites_idx],idx] =\
        #            (catalog[indices_of_prim_haloprop_bin[satellites_idx]]['halo_'+coord]
        #             -catalog[indices_of_prim_haloprop_bin[satellites_idx]]['halo_'+coord+'_host_halo']
        #            +catalog[rand_host_idxs]['halo_'+coord] + hc)%Lbox
        shuffled_pos[indices_of_prim_haloprop_bin[satellites_idx],idx] =\
            (catalog[rand_host_idxs]['halo_'+coord] + hc)%Lbox
            
    shuffled_upids[indices_of_prim_haloprop_bin[centrals_idx]] = -1
    shuffled_upids[indices_of_prim_haloprop_bin[satellites_idx]] = \
                catalog[rand_host_idxs]['halo_id']
        
    shuffled_host_mvir[indices_of_prim_haloprop_bin[centrals_idx]] = \
                catalog[indices_of_prim_haloprop_bin[centrals_idx]]['halo_mvir']
    shuffled_host_mvir[indices_of_prim_haloprop_bin[satellites_idx]] = \
                catalog[rand_host_idxs]['halo_mvir_host_halo']

for i, idx in enumerate(indices_of_prim_haloprop_bin[satellites_idx]): try: if shuffled_upids[idx] == catalog[idx]['halo_upid']: assert all([abs(shuffled_pos[idx, 0] - catalog[idx]['halo_x_host_halo']) < 1.0, abs(shuffled_pos[idx, 1] - catalog[idx]['halo_y_host_halo']) < 1.0, abs(shuffled_pos[idx, 2] - catalog[idx]['halo_z_host_halo']) < 1.0 ]) else: assert all([abs(shuffled_pos[idx, 0] - catalog[idx]['halo_x_host_halo']) > 1.0, abs(shuffled_pos[idx, 1] - catalog[idx]['halo_y_host_halo']) > 1.0, abs(shuffled_pos[idx, 2] - catalog[idx]['halo_z_host_halo'])>1.0]) except AssertionError: print i,idx print shuffled_upids[idx],catalog[idx]['halo_upid'] print shuffled_pos[idx] print catalog[idx]['halo_x_host_halo'], catalog[idx]['halo_y_host_halo'],catalog[idx]['halo_z_host_halo'] print '*'

from itertools import cycle colors = cycle(sns.color_palette()) for color, halo_id in zip(colors, np.unique(shuffled_upids[indices_of_prim_haloprop_bin[satellites_idx]])): idxs_of_sat_members = np.where(shuffled_upids[indices_of_prim_haloprop_bin]==halo_id)[0] idx_of_host = np.where(catalog[indices_of_prim_haloprop_bin]['halo_id'] == halo_id) plt.scatter(shuffled_pos[indices_of_prim_haloprop_bin[idxs_of_sat_members],0],\ shuffled_pos[indices_of_prim_haloprop_bin[idxs_of_sat_members],1], color = color, alpha = 0.8) plt.scatter(shuffled_pos[indices_of_prim_haloprop_bin[idx_of_host],0],\ shuffled_pos[indices_of_prim_haloprop_bin[idx_of_host],1], color = color, marker = 'v', s = 100, alpha = 1.0) plt.show()



In [13]:

    
catalog['halo_shuffled_vpeak_mag'] = shuffled_mags[:,0]
catalog['halo_shuffled_vvir_mag'] = shuffled_mags[:,1]
catalog['halo_shuffled_alpha_05_mag'] = shuffled_mags[:,2]
catalog['halo_shuffled_x'] = shuffled_pos[:,0]
catalog['halo_shuffled_y'] = shuffled_pos[:,1]
catalog['halo_shuffled_z'] = shuffled_pos[:,2]
catalog['halo_shuffled_upid']=shuffled_upids[:]
catalog['halo_shuffled_host_mvir'] = shuffled_host_mvir[:]



In [14]:

    
catalog.write('abmatched_halos.hdf5', format = 'hdf5', path = './abmatched_halos.hdf5', overwrite=True)



In [15]:

    
%%bash
ls -lt *.hdf5









    



-rw-r--r-- 1 swmclau2 des 557130076 Sep 27 11:48 abmatched_halos.hdf5



In [ ]:



In [ ]:



In [ ]: